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A NEW TIME-SPACE ACCURATE SCHEME FOR HYPERBOLIC PROBLEMS Is 

QUASI-EXPLICIT CASE* 

DAVID SIDILKOVERt 


Abstract. This paper presents a new discretization scheme for hyperbolic systems of conservation 
laws. It satisfies the TVD property and relies on the new high- resolution mechanism which is compatible 
with the genuinely multidimensional approach proposed recently. This work can be regarded as a first step 
towards extending the genuinely multidimensional approach to unsteady problems. Discontinuity capturing 
capabilities and accuracy of the scheme are verified by a set of numerical tests. 

Key words, unsteady problems, high- resolution schemes, TVD property. 

Subject classification. Applied and Numerical Mathematics. 

1. Introduction. A general approach towards the construction of genuinely multidimensional high- 
resolution schemes for hyperbolic systems was proposed recently ([16], [15], [17]). The main advantage of this 
approach is that it results in the discrete schemes that have some desirable properties which allow to design 
very efficient steady solvers (see [13] for a summary). However, these discretizations are suitable exclusively 
for steady problems since they are second order accurate at the steady-state only. 

Extending of the genuinely multidimensional approach to the transient problems may seem at the first 
sight straightforward. Noting that the unsteady advection problem in one dimensions is very similar to the 
steady advection problem in two dimensions, one can apply, say, the “fluctuation-splitting” (or “residual 
distribution”) steady 2D advection scheme for triangular meshes (proposed in [3] and given an algebraic 
formulation using limiters in [19]) to a one dimensional unsteady problem as well. However, a closer look 
reveals that the result of such a straightforward application may be that the “present” depends on the 
“future”. One could leave the causality concerns to the philosophers and try to iterate back and forth in 
space-time. However, the memory requirements of this approach will become absurdly demanding in more 
than one spatial dimension, since all the time instances need to be stored. 

This paper presents a new discrete scheme for the hyperbolic conservation laws and systems of conser- 
vation laws. It is one-step one-stage and satisfies the TVD property. The key feature of the new scheme is 
that it relies on the non-standard high-resolution mechanism which is compatible with the genuinely mul- 
tidimensional approach proposed recently and will allow to extend it to unsteady problems. The scheme is 
formulated using the framework introduced by Sweby [21], which makes it plain what arc the differences 
and similarities between the proposed and the existing discretizations. The new scheme has links to the 2D 
steady advection control volume scheme constructed in [14], [18], As the first stage we address only the case 
of Courant number less than one. The constructed scheme is not explicit even in this case. However, the 
treatment of the discrete equations can be made sufficiently similar to explicit (see §3). That is why we call 
the approach “quasi-explicit” . We verify the accuracy and the discontinuity capturing properties of the new 
scheme and demonstrate that the quality of the solutions obtained using it is at least as good as when using 
the standard schemes. 

*This work was supported by the National Aeronautics and Space Administration under NASA Contract Nos. NAS 1-1 9480 
and NAS1-97046 while the author was in residence at the Institute for Computer Applications in Science and Engineering 
(ICASE), NASA Langley Research Center, Hampton, VA 23681-2199. 

tlCASE, Mail Stop 403, NASA Langley Research Center, Hampton, VA 23681 
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The goal of this research direction is to construct very efficient implicit solvers for the unsteady com- 
pressible flow equations. The unsteady genuinely- multidimensional schemes are expected to allow this since 
they should maintain some fundamental properties of their steady ancestors: a good measure h-ellipticity of 
the implicit discretizations and the possibility to distinguish between different co- factors of the systems of 
equations. Therefore, the next step will be to consider the implicit case (Courant number larger than one) 
and to design an efficient implicit solver. This study will be reported elsewhere. 

The paper is organized as follows: §2 contains the description of the scheme (in control volume formu- 
lation) and the discussion concerning some basic properties. A possible solution procedure for the discrete 
equations is presented in §3. Generalization of the scheme to the Euler system is given in §4. §5 reports 
the numerical tests conducted with the constructed scheme. Some conclusions and discussion regarding the 
future research directions are presented in §6. The residual-distribution (fluctuation- splitting) formulation 
of the scheme is given in Appendix A. Also a brief discussion on more general conditions for a scheme to be 
TVD is given in Appendix B. 

2. The new scheme and its basic properties. Consider first an initial value problem for a scalar 
conservation law 

u t + f(u) x = 0, t > 0, x € K, 

u(x y 0) = Uo(#). 

Following Sweby [21] we shall use the following simplified notation (unless specified otherwise) for the data 
on the new and old time levels respectively 

(2.2) u k = u n k + \ u k = u n k . 


A general conservative discretization of (2.1) can be given as follows 


(2.3) 


u k = u k -X{h k+ L -h k _i) 


where h k + 1 is the consistent numerical flux 


(2.4) 


I — — mi ' ’ * j m+1 > ^ ) 

h(u, • • • , u) = f(u) 


and A is the mesh ratio 
(2.5) 


A = 


At 

Ax 


A crucial property needed for convergence proof of the difference scheme (2.3) is the so-called Total Variation 
Diminishing (TVD) property [4] 


( 2 . 6 ) 


TV{u n+1 ) < TV(u n ), 


where, abolishing temporarily the notation (2.2), we define the Total Variation an time level n as 
(2.7) 7 V(u") = £K +1 -<| 

n 

If the scheme (2.3) can be written in the form 


( 2 . 8 ) 


u k = u k — C fc _i Au fc _i +\iAu h ., 


1 I 
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then a sufficient condition for the scheme to be TVD (Harten’s Lemma, see [4]) is that the (data-dependent) 
coefficients satisfy the inequalities 

( 2 - 9 ) 0<C fc+ >, 0 <D k+ ,, 

and 


( 2 . 10 ) 


9 ^ C k+ 1 + D k+ i < 1. 


2.1. First order upwind scheme. Consider first for the purpose of illustration the first order upwind 
scheme. 

Define the local propagation speed 


ttk+i Au fc+ . 


( 2 . 11 ) 

where 

( 2 - 12 ) A/ fc+ i — /(ufc+i) - f{u k ), Au k+ i = u k +i - u k . 

A numerical flux defining a first order upwind scheme can be written as follows 


(2.13) 
Denote 

(2.14) 
and 

(2.15) 


fyfc+± — 2 l(f( U k) + f( u k + l)) l a *:4-i |Attfc + l]. 

a fc + i = |( a fc+£ + l a fc+i|) 
a fc+$ = l( a *+J “ l a fc+i|) 

V U\ = H + i 

vT i = Aa7 


'fc+i “ w fc+i' 
The local CFL number then can be written as follows 


( 2 ‘ 16 ) u *+h ~ Aa *+5 ~ v t+\ + u k+\- 

It is easy to see that using (2. 14), (2. 15) the scheme defined by (2.13) will read as 

(2.17) 

Denoting 


u = Wfc - *£_ i Au fc+ 1 


(2.18) 




D t , i = -j/, 




makes it obvious that inequality (2.9) is satisfied, and that 
(2.19) C l 

is a CFL condition. Therefore, the first order upwind scheme is TVD 


'k+i + Dk+% — v + 
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Fig. 1 . Stencil of the Lax-Wendroff scheme. 


2.2. The high-resolution scheme: constant coefficient equation. In order to illustrate the con- 
struction of a high-resolution scheme we shall consider first a linear constant coefficient equation 


( 2 . 20 ) 


u t 4- au x = 0, a > 0. 


The first order upwind scheme approximating (2.20) can be written as follows 


( 2 . 21 ) 


where 


u k = Uk — vAu k _ i 


is the CFL number. The second order Lax-Wendroff [6] scheme can be given as follows (the corresponding 
stencil is depicted on Fig.l) 


(2.23) 


u k - - A_[-(l - i/)vAu k+ i] 


or as a first order upwind scheme (2.21) whose flux (2.13) is augmented by the following anti-diffusive flux 


(2.24) 


— 2^(1 - ^Au k+k 


2.2.1. Some existing schemes. The Lax-Wendroff (as well as any second order linear scheme) is 
known not to satisfy the TVD property. A high- resolution scheme can be defined by adding a “limited” 
anti- diffusive flux to the first order upwind one. The resulting scheme reads as follows 


(2.25) 


u k - vAu k _i - A-{xp(q k )-(l - is)i/Au k+ i], 


where q k is standardly defined as a ratio of two subsequent finite differences of the numerical solution on the 
time- level n 


Noting that the following identity holds 


(2.27) 


1 p{q k ) Au k+ 1 = 


it can be concluded (see [21]) that high-resolution scheme (2.25) can be viewed as a nonlinear hybrid of Lax- 
Wendroff and the Warming-Beam [22] schemes. The latter is given by the following (the stencil is depicted 
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Fig. 3. TVD region. 


at Fig.2) 

(2.28) u k = Uk — i/A u k _i — A_[^(l — v)vAu k _ 1 ] 

2 2 2 

It was also shown in [21] that the high-resolution scheme (2.25) is TVD provided the limiter function ^ 
satisfies the following requirements 

(2.29) 0 < ip(q) <2, 0 < < 2 

9 

The common practice is also to require 

(2.30) = 0, q < 0. 

Hence (see [21]) if the limiter function lies within the “shaded” area shown on Fig.3, scheme (2.25) satisfies 
TVD property. Requiring the scheme also to be second order accurate wherever possible implies that 'ip(q) 
has to be a Lipschitz-continuous function and 

(2.31) V’(l) = 1- 

It is also desirable (see [21]) that the resulting non-linear high-resolution scheme is an internal average of 
the Lax-Wendroff and Warming-Beam scheme. Therefore, we obtain the second order subset of the TVD 
region as depicted on Fig. 4. 


5 




Fig. 4. Second order TVD region. 


2.2.2. The new high-resolution scheme. Define the argument of the limiter function in the following 


way 

(2.32) 

where 


A t Ufc 


(2.33) 


A t Ufc = u k —Uk 


is a finite difference of the numerical solution values at times n-f 1 and n at point h. In other words, the 
limiter function in the new scheme relies on the ratio of the finite differences in time and space (each taken 
with an appropriate coefficient). 

The identity 

(2.34) %p{q k )i^Au k+ k = tUk 

allows to rewrite scheme (2. 25), (2. 32) in the following form 

(2.35) u k = u k - vAu k -i - 1(1 - v)[- r ^-A t u k + A t u fc -i], 

suggesting that the scheme (2. 25), (2. 32) (whose stencil is represented by Fig. 6) can be viewed as a nonlinear 
hybrid of the Lax-Wendroff and the “box” schemes. The latter is given by the following 

(2.36) u k = u k - vAu h _i - 1(1 - i/)[-A t u fc + A t u fc _!] 


and its stencil is depicted on Fig. 5. 

Using (2.34), we can rewrite (2. 25), (2. 32) as follows 


(2.37) 


u k - uA u k _^ - 1(1 - v)[-~^A t u k - V’(<i , Jt-i)^A« fc _ 1 ], 


Qk 


or 

(2.38) 


u k = u k - 


1 ~ ~ ^(gfc-l) . 

-*#(»)/» U ' 


2 


;l I 
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Fig. 5. Stencil of the “Box” scheme. 



Fig. 6. Stencil of the new high-resolution scheme. 


The latter can be written in the form (2.8) with 


(2.39) 

and 


c i = 1 - - ^(gfc-i) 

1 ~ l(l 


(2.40) Dk+i = 0 
It can be easily verified that the inequality (assuming v < 1) 

(2.41) 0 < Cfc_i < 1 
is equivalent to the following 

(2.42) rp(q) < , and ~ vil>{q) < 2, 

which, in its turn, is satisfied if (2.29) holds. 

Remark 2.1. The idea of using the %ox ” scheme as one of the building blocks when constructing a 
discretization for hyperbolic problem can be found in [9]. It is employed for the case of CFL > l, while the 
Lax- Wendroff scheme is used when CFL < 1. 

2.3. General case. The new high-resolution scheme for general nonlinear problem (2.1) can be defined 
by adding a limited anti- diffusive flux to that of the first order upwind scheme (2.17) (as it was done previously 
for the linear constant coefficient case) 

( 2 - 43 ) h k+\ = K+i + - v k+ i) )«fc+i + V’fafc+iK+i) Au *+i> 
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where 


(2.44) 


+ _ 

Qk Xa t+i Au k+ i 


and 

(2.45) 


Qk- i-i — 




Aa, . ! Au 




Using the identities 

(2.46) 
and 

(2.47) 


,, + x + A _ V’(9fc)A f w fe 

WK+iA^+i = -f t- 


V’tefc+iK+^Aufc+i 


V'fafc+i) A t Ufc +1 


9fc+i A 

the alternative form of the numerical flux (2.43) can be obtained 


(2.48) 

Define 

(2.49) 
and 

(2.50) 
Denoting 

(2.51) 
and 

(2.52) 


K+k = 


,=K +l - hi - ^ +l) ( + Mai A ‘ u * +1 ) . 
2 fc+5 2 2 V 4 A 9fc+i A ) 




^ + i = ^(9fc + iK + x(l - " fc+J )> A* fc+ | = - V l + J- 


iC 1 

. " 2 Z 2 




2 i-Vl-l-Vlr+l 


n _ 

^ k ^r 4 1 - 4. 

-#**+* 


the scheme (2.3) with numerical fluxes defined by (2.43) (or (2.48)) can be put into form (2.8). 

Finally, we can state the following result 

Theorem 2.2. Scheme (2.3) with high resolution flux (2.43) that incorporates a limiter function satis- 
fying (2.29) is TVD. 

Proof. Note that if v k _\ > 0, then fi k _i = 0. Also, if i/+ +i > 0, then j. i* +i = 0. It follows from here 

that 


(2.53) C k ^> 0, £> fc+ i> 0, 

provided the limiter function tp satisfies inequality (2.29). 


8 



Inequality (2.29) gives the following 
(2.54) < 1, D k+ , < 1, 

However, since either or l/ k+ i vanishes, either C fc+ i or D k +± w ill vanish as well. Therefore, 

( 2 - 55 ) C kH +D k+l < 1, 

which together with (2.53) completes the proof. □ 

Remark 2.3. Notice that a purely upwind decomposition (2-44)~(^4 5 ) w necessary here for the con- 
struction , as the ratio between the numerator and denominator must be close to 1 in smooth monotone regions 
in order to keep accuracy. In particular , a Lax- Friedrichs type splitting cannot be used here. 

3. Solution procedure. Consider for simplicity scheme (2.25) with (2.32), devised for linear problem 
(2.20). This scheme is not explicit the stencil involves more than one point on the new time level n- hi 
(Fig. 6). However, there is no need for a global solver in this case. One can use some simple local procedures 
for solving procedures for solving the resulting discrete equations. 

Space marching . 

An obvious possibility is a space marching, i.e. solve the discrete equation at point k after it has been 
solved at point k — 1. However, considering solving the nonlinear systems of equations and/or extensions to 
multi- dimensions, it appears that the latter method may be not always feasible/desirable. 

Predictor- Corrector. 

Consider an initial-boundary value problem for (2.20) discretized on a grid k = 0 Assume, that 
Ufc, the solution on the time level n, is a monotonic function. 

Consider an arbitrary discrete function as an initial guess for a solution at the time level n + 1 . 
Performing one sweep of Jacobi relaxation using scheme (2.25) with (2.32) will result in u fc *, i = 0, ..., K. 
We can state the following 

Lemma 3.1. If the discrete function u k , k = 0, ...,K is monotonic , then the discrete function u k * is 
monotonic as well . 

Proof Putting scheme (2.25), (2.32) into form (2.8) with C k _±,D k+ 1 defined by (2. 39), (2. 41) and (2.39) 
respectively, we can conclude that 

(3.1) min(ufc_i,ufc) < u k * < max^-i, u k ), 

which concludes the proof. □ 

This suggests that the following predictor- corrector procedure can be used for solving the discrete equa- 
tions resulting from the limited scheme: 

Predictor. Lax-Wendroff scheme to produce a grid-function u k \ k = 0,...,Ff, which is a second order 
accurate solution to (2.1), but may exhibit an oscillatory behavior near discontinuities; 

Corrector. Perform a Jacobi relaxation to update values u k t to obtain a new grid-function u k , k = 0, ..., K 
that is non-oscillatory (due to Lemma 3.1) and second order accurate solution to (2.1). 

Such a solution procedure is essentially a point-implicit solver. Note that the resulting discrete function may 
not in general satisfy the discrete equations corresponding either to the limited scheme or to the Lax-Wendroff 
scheme. 

This should not cause any problems in the smooth regions since the discrepancy between the solutions 
to two different second order schemes is expected to second order small. 
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The difficulty though may arise in the vicinity of discontinuities. Different schemes may result in different 
discontinuity profiles. Therefore, if the residuals of one scheme are measured on the solution corresponding 
to another scheme, they may be large. If no special care is taken, this may have a cumulative effect over 
several time-steps and lead to wrong shock speeds. 

Our experience though demonstrates, that it is sufficient to perform a few (instead of one) updates of 
the solution in the vicinity of a discontinuity using the “limited” scheme to make sure that the residuals of 
this scheme are very small. We shall elaborate on this issue in the future publication which will address also 
the implicit case. 

4. Hyperbolic systems. We shall briefly sketch out an extension of the new scheme to the hyperbolic 
systems of conservation laws, in particular, to the Euler equations of gas dynamics 

(4.1) u t + f(u) x =0 

Again, as in the scalar case, the discretization can be put into the following general form 

(4.2) u fe = u fc - X(h k+ 1 - h fc _ i ) 


(4.3) u* = u£ +1 , u fc = u£ 

For the case of the Euler equations we can take the Roe-averaged Jacobian matrix of f (see [11]) 

(4.4) A k+ i (ufc +1 - Ufc) = f(ufc +1 ) - f(ufc) 
and define the first order upwind flux as follows 

hfc + i = h(f(u fc ) + f(ufc +1 )) - |ifc + i|(ufc +1 - Ufc)] 


(4.5) 

Here 

(4.6) |ifc + i| = T|A|T-\ 

T is the matrix of right eigenvectors of A, and 

(4-7) A = T- l A k+l T, 

is a diagonal matrix 

(4.8) A = diag{ai,...,a n } 
and its absolute value is defined as 

(4.9) |A| = diag{|oi|, ..., |a„|}. 

We define the flux corresponding to the new high-resolution scheme as 

(4-10) 

where the anti-diffusive flux is given by 

(4.11) h AD = (h l ,...,h n ) T 


h fc+ i = K + i + \T~ l K+ h 
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and 

(4.12) 
with 

(4.13) 


( h i)k+ i = C 1 ~ ("«)*+$) (V'((9i)fc)(ai)fc + i + V’((9i) fc+ i)(a«) fc+ i) A(wi) k+ i 


(*) 


k-\- 


a iAt 
5 = ~Ax t 


(4.14) 


<«>: - — ^k. 


A(a i )+ +i A(w; i ) fe+ i ’ 


(4.15) 


(?')*+! — 




and (i /»)+ + 1 and (i^)” + A are the “positive” and “negative” parts of (i^) fc+ i respectively. 

Remark 4.1. IVbte i/iat the formulated scheme (as well as any other scheme based on the pure upwind 
decomposition) may produce non-physical solutions violating the entropy condition. Some ways to fix this 
problem have been described in the literature (see, for instance , [5]). Therefore , we do not include a discussion 
on this matter here . 


5. Numerical results. In this section we present some numerical tests of the new scheme using a few 
standard test problems. 

5.1. Scalar case. 


5.1.1. Linear advection equation. The purpose of this test is to verify the second order accuracy of 
the method. Therefore, we use here a simple problem with the known exact solution. 

Consider the initial- boundary value problem for the following equation 


(5.1) 


u t + .5 u x = 0. 


on the domain [— 1, 1] with the initial condition given by 

(5.2) u(x, 0) = — .5sin7rx + .5 


and the boundary condition (at point x — — 1) 

(5.3) u(— l,t) = ,5sin(l 4* .5£) + .5. 

It is obvious that the exact solution to this problem is given by 

(5.4) u = — .5sin7r(x — .5 1) + .5. 

The L\ error norm of the solution error at time t = 1.0 for different levels of resolution is presented in Table 
5.1. The Courant number in all the tests is .5. The error reduction due to the doubling of the number of 
the grid-points is roughly by factor 4, which verifies that the scheme is second order accurate for smooth 
solutions. 
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Table 5.1 

Linear advection equation: solution error L\ norm behavior with the mesh refinement. 


Number of grid-points 

Solution error {L\ norm) 

Numerical order of accuracy 

25 

0.0075311723 


50 

0.0018290140 

2.042 

100 

0.0004250877 

2.106 

200 

0.0000971690 

2.129 


Table 5.2 

Burgers equation: solution error L\ norm behavior with the mesh refinement . 


Number of grid-points 

Solution error ( L\ norm) 

Numerical order of accuracy 

50 

0.0020134142 


100 

0.0005068036 

1.990 

200 

0.0001261164 

2.007 


5 . 1 . 2 . Burgers equation. Our next purpose is to examine the accuracy of the new scheme as well as 
its capability to resolve shocks, representing them by a sharp oscillation-free profile. The problem considered 
here is the following nonlinear conservation law 

(5.5) u t + .5(u 2 ) x = 0. 

The following expression (raised sine) serves as the initial condition 

(5.6) u{x , 0) = — sin(7nr) + 0.5 

First we demonstrate the accuracy of the scheme for time t = .25 by comparing the solution error norm for 
different levels of resolution (see Table5.2). In this case the solution is still smooth, since the shock did not 
start to develop yet. The ability of the scheme to resolve shocks is demonstrated on the same test problem. 
The numerical solution to this problem corresponding to time t = 2 is presented on Fig. 7. We can see that 
the shock is represented by a sharp layer and that there are no over- and undershoots on its sides. 

5 . 2 . Euler equations. In this section we present some numerical tests of the new scheme using a few 
standard test problems concerning the Euler equations of gas dynamics. The scheme used is as described in 
§4 with the Van Albada limiter. Note, that the Van Albada limiter does not satisfy the inequalities (2.29) 
(see Fig. 12), which are a sufficient condition for a (scalar) scheme to be TVD. A more general sufficient 
condition, that is satisfied by the Van Albada limiter, is presented in Appendix B. The reasons for using the 
Van Albada limiter is discussed in Appendix B as well. 

All the test-cases considered here use the domain V — {x : x G [—1,1]}. 

Lax’s problem. The initial data for this problem are two constant states, one to the left of the origin 
(— 1 < x < 0): pi = .445, qi — .698, Pi = 3.528, another to the right of the origin (0 < x < 1): p r = .5 ,q T — 
0 ,P r = .571. Here p stands for density, q for velocity and P for pressure. 

The numerical results for this case are presented on Fig. 8. The computational grid is 200 points. Solid 
lines represent the exact solution to the problem. 

Sod’s problem . Similarly to the previous problem, the initial data are given for the two constant states: 
pi = l,qi = 0, Pi = 1 and p r = .125 ,q r = 0,P r = .1. 

The grid consists of 200 points and the solid line represents the exact solution. 
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Fig. 7. Burgers equation. A numerical solution obtained on a grid of lOOpts. at time t ~ 2. The initial condition is given 
by (5.6). 

Shock interacting with a disturbance in density (the testcase by Shu & Osher [12]). The initial data for 
this problem is pi - 3.857143, qi = 2.629367, P t = 10.33333 for -1 < x < -.8 and p r = 1 + .2sin(57r x),q r = 
0,P r = 1 for -.8 < x < 1. 

The interaction of the moving shock with a sinusoidal disturbance in density generates higher frequencies 
in the post- shock area. Here we test a capability of the new scheme to resolve these oscillations. The 
numerical experiments use three different levels of resolution (see Fig.l0a,b,c). The solid line represents a 
numerical solution obtained on a very fine grid (1600 points). The snapshots correspond to time t = .36. 

Colella- Woodward problem . The initial data for this problem is given by three constant states: pi = 
1 ,qi = 0 ,P t = 1000 for -1 < x < -.8, p m = l,g m = 0, P m = .01 for -.8 < x < .8 and p r = l,g r = 0, P r = 
100 for .8 < x < 1. Results are shown on Fig. 11 (time t = .75). 

6. Conclusions. A new discretization scheme for hyperbolic problems was presented. It can be viewed 
as high- resolution generalization of the Lax-Wendroff scheme. In this paper we consider only the case of 
Courant number less than unity. Even in this case the resulting scheme is implicit. However, we demonstrate 
that the solution procedure in this case does not have to involve a implicit global solver, but can rely on 
some simple local techniques. That is why we call this case quasi- explicit 

The methodology incorporated into the new high- resolution mechanism should allow to extend the 
previously proposed genuinely multidimensional approach (see [16, 15, 17]) to unsteady problems. This is a 
subject of the future study. 

Another very promising direction is to generalize the scheme for the case of Courant number larger than 
unity. The procedure of solving the discrete equations in this case may require a global implicit solver. 

The properties of the genuinely multidimensional (steady) schemes summarized in [13], however, indicate 
that extremely efficient implicit solvers can be constructed for unsteady generalizations of such schemes as 
well. This is another major direction for the future research. 

It is also interesting to note that the schemes proposed by Colella [1, 2], LeVeque [7], [8] and Radvogin 
[10] can all be regarded as extensions to multi-dimensions of scheme (2.25) with the standard limiters with 
the argument given by (2.26). The genuinely multidimensional scheme (in the spirit of [15] when extended 
to unsteady problems can be viewed as extension of the scheme presented in this paper, namely, (2.25) with 
the limiters relying on the ratio of differences in space and time (2.32). 
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Fig. 8. Lax’s problem. A numerical solution obtained on a grid of 200pts. is represented by circles . The solid line 
corresponds to the exact solution. 

Consider a transient problem that converges to a steady-state. It is interesting to note that if the Lax- 
Wendroff scheme (or one of its standard generalizations) is applied to solve such a problem, the computed 
steady-state solution will depend on the time-step (since the artificial dissipation of such schemes scales with 
the time-step). The scheme proposed in this paper, however (even though it is an extension of the Lax- 
Wendroff scheme), will produce a steady-state solution that is independent of the time-step. This is due to 
the fact that if the temporal changes of the solution become very small and eventually vanish, the limiter will 
cut off all the spatial terms in the discretization that scale with the time-step. Thus, the scheme essentially 
reduces to the first order upwind. The generalization of the scheme presented here to multi- dimensions is 
expected to maintain this property to a certain extent as well. 
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Appendix A. Residual-distribution (fluctuation-splitting) formulation. 

The approach for constructing genuinely multidimensional upwind schemes for the Euler equations was 
formulated in [16], [15] in the residual-distribution context. Here we shall reformulate the scheme constructed 
in §2.3 in the same context as well. 

The philosophy of this approach, proposed in [3] for scalar advection is that the discrete equation to 
be solved at each grid point is constructed from portions of the residuals (fluctuations) of the equation(s) 
evaluated on the grid elements (segments in one-dimensional case) adjacent to this grid point. The question 
of constructing a discrete scheme reduces to the question of defining a rule of distributing residuals to the 
nodes from a grid element. 

A.l. Scalar case. Consider a scalar conservation law (2.1). Its residual on the segment [k y k + 1] 
(multiplied by the mesh- size) is given by 

(A.l) r = A/ fe+ i = a fc+ 1 A« fc+ i 

The first order upwind scheme at the points k and k - hi can given by the following residual distribution 
formulas 

A t Uk+i = (r + r) + C.O.L. 

[ ’ A t u k = \\{r-f) + C.O.L. 

where C.O.L . stands for Contributions from Other Elements and the term f stands in this case for the 
artificial dissipation and is defined as follows 


(A.3) r = sign(a fc+ 1 )r = |a fc+ i|Aw fc+ i 

The higher-resolution scheme identical to that presented in §2.3 can be obtained by defining 


(A.4) 


f = sign(a fc+ i) [r + (1 - v k+ ±) (rp(q^)a+ + ^ + V’fafc+iK+i) Au fe+ i 


where the quantities are as defined in §2.3. 

A. 2. Euler system. The residual of the Euler system on the segment [A:, A: -f 1] can be computed 
according to the following 


( A.5) r = Af fc+ 1 = A fc+ 1 Au fc+ i 

By analogy to the scalar case, the first order upwind scheme is defined by the residual distribution formulas 


(A.6) 

Au fc+1 — |A(r + r) + C.O.L. 
Au fc = |A(r - f) + C.O.L. 

with 


(A.7) 

f = sign(A fc+ 1 )r = | A k+ 1 \ Au i+ ^ 


The high-resolution scheme identical to one constructed in §4 can be obtained by the following definition 
(A. 8) r = sign(i fc+ i ) [r + T -1 h£°jj . 

Again, the notation used in (A. 7) and (A. 8) was established in §4. 
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Fig. 12. Van Albada limiter . 


Appendix B. Van Albada limiter and TVD property. Implicit solvers for steady-state or unsteady 
computations based on the high- resolution schemes are usually unable to drive residuals to the machine zero. 
The reason for this is the numerical fluxes that may be discontinuous due to the limiters. Limiters are based 
on the ratio of two quantities. A limiter function can be looked at as a function of two arguments: the 
numerator and denominator of this ratio. The limiter function that satisfies the “standard” TVD condition 
(see (2.29 or Fig. 3) is a discontinuous function of its second argument (denominator of the ratio) at zero. 

It is not clear whether or not the inability to obtain convergence of the residuals to the machine zero is 
a serious problem, since the “limit cycle” normally occurs when convergence below the level of truncation 
error is achieved. However, this difficulty can be easily removed by using a limiter function which satisfies 
the following requirement 


(B.l) 


t^l-hOG [ — OC ’ 


For instance, the Van Albada limiter, which is given by 


(B.2) 


ip{q) = 


q 2 + q 

q 2 + 1 


(see also Fig. 12), satisfies (B.l). The problem, however, seems to be that it does not satisfy the “standard” 
TVD condition (2.29). 

Note that (2.29) is a sufficient condition for a high- resolution scheme to be TVD and, in fact, is too 
restrictive. A more general condition on the limiter in order to ensure the TVD property of the scheme is 
presented, for instance, in [20]. 

Here we shall just state that the following generalization of (2.29) 


(B.3) —a < ip(q) <2 - a, -a < < 2 - a 

q 

with a > 0, also implies the TVD property. The latter is quite clear from (2.42) in §2.2.2. A necessary 
condition for the second order accuracy is a < 1. 

It is obvious that the Van Albada limiter satisfies (B.3) with, for instance, a = .5. 
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